Signals are calculated using the Shockley-Ramo theorem. 
The current \(i\left(t\right)\) induced by a particle with charge 
\(q\) at a position \(\mathbf{r}\) moving at a velocity \(\mathbf{v}\)
is given by
\begin{equation}\label{Eqn:RamoShockleyWeightingField}
  i\left(t\right) = -q \mathbf{v} \cdot \mathbf{E}_{w}\left(\mathbf{r}\right),
\end{equation}
where \(\mathbf{E}_{w}\) is the so-called weighting field for the 
electrode to be read out and the charge induced by particle moving from 
$\mathbf{r}_{1}$ to $\mathbf{r}_{2}$ is given by
\begin{equation}\label{Eqn:RamoShockleyWeightingPotential}
  \int\limits_{t_{1}}^{t_{2}}i\left(t\right)\text{d}t = q\left[\phi_{w}\left(\mathbf{r}_{2}\right) - \phi_{w}\left(\mathbf{r}_{1}\right)\right],
\end{equation} 
where $\phi_{w}$ is the weighting potential.

The basic steps for calculating the current induced 
by the drift of electrons and ions/holes are:
\begin{enumerate}
  \item
  Prepare the weighting field and/or weighting potential 
  for the electrode to be read out. 
  This step depends on the field calculation technique 
  (\textit{i.\,e.} the type of \texttt{Component}) that is used 
  (see Chapter~\ref{Chap:Components}). 
  \item
  Tell the \texttt{Sensor} that you want to use this 
  weighting field/potential for the signal calculation. 
\begin{lstlisting}
void Sensor::AddElectrode(ComponentBase* cmp, std::string label);
\end{lstlisting}
  where \texttt{cmp} is a pointer to the \texttt{Component} 
  which calculates the weighting field/potential, and \texttt{label} 
  (in our example \texttt{"readout"}) is the name 
  you have assigned to the weighting field/potential in the previous step.
  \item
  Setup the binning for the signal calculation.
\begin{lstlisting}
void Sensor::SetTimeWindow(const double tmin, const double tstep, const int nbins);
\end{lstlisting}
  The first parameter in this function is the lower time limit (in ns), 
  the second one is the bin width (in ns), and the last one 
  is the number of time bins.
  \item
  Switch on signal calculation in the transport classes using 
\begin{lstlisting}
void AvalancheMicroscopic::EnableSignalCalculation();
void AvalancheMC::EnableSignalCalculation();
void DriftLineRKF::EnableSignalCalculation();
\end{lstlisting}
  The \texttt{Sensor} then records and accumulates the signals of all 
  avalanches and drift lines which are simulated.
  \item
  The calculated signal can be retrieved using 
\begin{lstlisting}
double Sensor::GetSignal(const std::string label, const int bin);
double Sensor::GetElectronSignal(const std::string label, const int bin);
double Sensor::GetIonSignal(const std::string label, const int bin); 
\end{lstlisting}
  The functions \texttt{GetElectronSignal} and 
  \texttt{GetIonSignal} return the signal induced by negative 
  and positive charges, respectively. \texttt{GetSignal} returns 
  the sum of both electron and hole signals.   
  \item
  After the signal of a given track is finished, call
  \begin{lstlisting}
void Sensor::ClearSignal();
  \end{lstlisting}
  to reset the signal to zero.
\end{enumerate}

For plotting the signal, the class \texttt{ViewSignal} can be used.
As an illustration of the above recipe consider the following example. 
\begin{lstlisting}
// Electrode label
const std::string label = "readout";
// Setup the weighting field.
// In this example we use a FEM field map.
ComponentAnsys123 fm;
...
fm.SetWeightingField("WPOT.lis", label);

Sensor sensor;
sensor.AddComponent(&fm);
sensor.AddElectrode(&fm, label);
// Setup the binning (0 to 100 ns in 100 steps).
const double tStart =   0.;
const double tStop  = 100.;
const int nSteps = 100;
const double tStep = (tStop - tStart) / nSteps;

AvalancheMicroscopic aval;
aval.SetSensor(&sensor);
aval.EnableSignalCalculation();
// Calculate some drift lines.
...
// Plot the induced current.
ViewSignal signalView;
signalView.SetSensor(&sensor);
signalView.PlotSignal(label);
\end{lstlisting}

The algorithms used for calculating the induced current are slightly different 
for \texttt{DriftLineRKF}, \texttt{AvalancheMC}, and \texttt{AvalancheMicroscopic}.

For drift lines calculated using the RKF method, 
the times $t_j$, coordinates $\mathbf{r}_j$ and drift velocities $\mathbf{v}_j$
at each point along the drift line are taken and the induced current 
\begin{equation*}
  i_{j} = -q \mathbf{E}_{w}\left(\mathbf{r}_{j}\right) \cdot \mathbf{v}_j
\end{equation*}
at these points is computed. In order to calculate the average current 
in each time bin, the array of $\left(t_{j}, i_{j}\right)$ is interpolated 
(linearly) and then integrated using Simpson's rule over $2n_\text{avg} + 1$ points. The parameter $n_\text{avg}$ defaults to 2 and can be set using
\begin{lstlisting}
void DriftLineRKF::SetSignalAveragingOrder(const unsigned int navg);
\end{lstlisting} 

For drift lines simulated using microscopic tracking 
(\texttt{AvalancheMicroscopic}), the signal between subsequent points 
along the drift path is assumed to be constant, and can be 
calculated using either the weighting potential or the weighting field. 
The method can be selected using 
\begin{lstlisting}
void AvalancheMicroscopic::UseWeightingPotential(const bool on = true);
\end{lstlisting}
By default, the weighting field is evaluated at the mid-point of a 
drift line segment. 
Using 
\begin{lstlisting}
void AvalancheMicroscopic::EnableWeightingFieldIntegration(const bool on = true);
\end{lstlisting} 
one can request 6-point Gaussian integration of the weighting field 
over a drift line segment.

The weighting potential method takes $\phi_{w}$ 
at the start and end of each step $j \rightarrow j + 1$ along the drift path
and calculates the average current $q\Delta\phi_{w}/\left(t_{j+1} - t_{j}\right)$.

For drift lines simulated using \texttt{AvalancheMC} one also has the choice 
between using the weighting potential or the weighting field for computing 
the induced current. The method to be used can be selected using
\begin{lstlisting}
void UseWeightingPotential(const bool on = true);
\end{lstlisting} 
The weighting potential method should be used if one wants to ensure that 
the integrated current is correct (equal to the collected charge).
\section{Readout electronics}

In order to model the signal-processing by the front-end electronics, the 
``raw signal'' -- \ie~the induced current -- 
can be convoluted with a so-called ``transfer function'' (often also referred 
to as delta response function). 
The transfer function to be applied can be set using
\begin{lstlisting}
void Sensor::SetTransferFunction(double (*f)(double t));
\end{lstlisting}
where \texttt{f} is a function provided by the user, or using
\begin{lstlisting}
void Sensor::SetTransferFunction(const std::vector<double>& times,
                                 const std::vector<double>& values);
\end{lstlisting}
in which case the transfer function will be calculated by 
interpolation of the values provided in the table.
A third option is to use a predefined expression, implemented in 
the helper class \texttt{Shaper}. Its constructor,
\begin{lstlisting}
Shaper(const unsigned int n, const double tau, const double g, std::string shaperType);
\end{lstlisting}
takes four arguments: $n$ is the order of the shaper, 
$\tau$ is the time constant, $g$ is the gain factor, 
and \texttt{shaperType} is either \texttt{"unipolar"} or \texttt{"bipolar"}.
In the first case (unipolar shaper), the transfer function is given by
\begin{equation*}
  f\left(t\right) = g \exp\left(n\right)\left(\frac{t}{t_{p}}\right)^{n}
                    \exp\left(-t / \tau\right), \qquad t_{p} = n\tau,
\end{equation*} 
while for a bipolar shaper the expression
\begin{equation*}
  f\left(t\right) = g \frac{\exp\left(r\right)}{\sqrt{n}} \left(n - \frac{t}{\tau}\right)
                    \left(\frac{t}{t_{p}}\right)^{n - 1} 
                    \exp\left(-t / \tau\right), \qquad t_{p} = r\tau, \qquad r = n - \sqrt{n}.
\end{equation*}
is used. The normalization of these expressions is chosen such that the  
value of the transfer function at the peaking time $t = t_{p}$ is unity.
In order to use a transfer function provided by a \texttt{Shaper} class,
one should call
\begin{lstlisting}
Sensor::SetTransferFunction(Shaper& shaper);
\end{lstlisting}

The presently stored signal can be convoluted with the 
transfer function (specified using one of the methods above) using 
\begin{lstlisting}
bool Sensor::ConvoluteSignal();
\end{lstlisting}

As an example, consider the following transfer function
\begin{equation*}
  f\left(t\right) = \frac{t}{\tau}\exp\left(1 - t/\tau\right), \qquad
  \tau = 25\,\text{ns},
\end{equation*}
\ie~a unipolar shaper with $n = 1$. The two code snippets 
below illustrate different methods for applying this transfer function 
to the induced signal. 
\begin{itemize}
\item 
\begin{lstlisting}
double transfer(double t) {
  constexpr double tau = 25.;
  return (t / tau) * exp(1 - t / tau);
}

int main(int argc, char* argv[]) {

  // Setup component, media, etc.
  // ...
  Sensor sensor;
  sensor.SetTransferFunction(transfer);
  // Calculate the induced current.
  // ...
  // Apply the transfer function.
  sensor.ConvoluteSignal();
  // ...
}
\end{lstlisting}

\item

\begin{lstlisting}
int main(int argc, char* argv[]) {

  // ...
  Shaper shaper(1, 25., 1., "unipolar");
  Sensor sensor;
  sensor.SetTransferFunction(shaper);
  // ...
  sensor.ConvoluteSignal();
  // ...
}
\end{lstlisting}
\end{itemize}

\subsection{Noise}
Prior to convoluting the induced current with a transfer function, 
one can add a random noise component to the signal using 
\begin{lstlisting}
void AddWhiteNoise(const double enc, const bool poisson = true, const double q0 = 1.);
\end{lstlisting}
\begin{description}
  \item[enc] desired equivalent noise charge (ENC) of the convoluted signal,
  \item[poisson] flag whether to sample the number of noise pulses 
    from a Poisson distribution, or to sample the noise charge in each 
    bin from a Gaussian distribution,
  \item[q0] amplitude of the noise delta pulses (in electrons). 
\end{description} 

The algorithm is based on the fact that white current noise is equivalent 
to a random sequence of delta current pulses with large frequency $\nu$ 
and with constant amplitude $q_0$. When processing this signal by an 
amplifier with transfer function $f\left(t\right)$, 
with a peak normalized to unity, the variance of the output signal becomes
\begin{equation*}
  \nu q_{0}^{2} \int f\left(t\right)^{2} \text{d}t.
\end{equation*}
We want this to be equivalent to the ENC$^2$, 
which defines 
\begin{equation*}
\nu = \frac{1}{q_{0}^{2}}\frac{\text{ENC}^2}{\int f\left(t\right)^2 \text{d}t}.
\end{equation*}

The total number of current delta pulses in a period of time $\Delta{t}$ 
is Poisson distributed with a mean $\nu\Delta{t}$ 
and a standard deviation $\sqrt{\nu\Delta{t}}$.
With the flag \texttt{poisson} set to \texttt{true}, a Poisson-distributed 
number of current pulses is added to the signal in the time window.

For large frequencies, the Poisson distribution becomes a 
Gaussian distribution, so the standard deviation of charge 
in a time bin $\Delta{t}$ is 
\begin{equation*}
  \sigma_{q} = q_{0} \sqrt{\nu\Delta{t}}.
\end{equation*}
With the flag \texttt{poisson} set to \texttt{false}, 
a Gaussian-distributed noise charge is added to each signal bin. 
